/*
 *	$Id: plotter.cpp 171 2009-06-21 22:46:53Z wabsie $
 */

#include "plotter.h" // Hvis andre #include's i denne fil end denne, så slet dem


AG_PixelFormat* plotter::format = AG_PixelFormatRGB(32, colorRed, colorGreen, colorBlue);

plotter::plotter(view par) :
	param(par),

	// Sørg for at vores surface ikke bliver mindre end 10x10 pixels
	w((par.plotWidth > minimumPixels)  ? par.plotWidth  : minimumPixels),
	h((par.plotHeigth > minimumPixels) ? par.plotHeigth : minimumPixels),

	// Sæt view-range:
	xmin(par.firstAxMin), xmax(par.firstAxMax),
	ymin(par.secondAxMin), ymax(par.secondAxMax),
	xrange(xmax - xmin),
	yrange(ymax - ymin),

	centre(w/2, h/2),					// Pixel-coords til pixmap'ets centrum
	scale(w/xrange, -h/yrange),
	offset(xmax+xmin, ymax+ymin),
	origo(centre - scale.mul(offset/2))	// Pixel-coords til origo
{
	AG_TextColor32(colorDefault);		// Sæt tekstfarven

	// Opret en surface vi kan tegne på
	this->plot = AG_SurfaceNew(w+1, h+1, format, AG_SWSURFACE);
    this->skala = AG_SurfaceNew(par.intensityWidth, par.intensityHeigth, format, AG_SWSURFACE);
	this->zbuffer = 0; // Z-buffer allokeres senere, hvis nødvendigt
}

plotter::~plotter()
{
	if (this->zbuffer != 0)
		delete this->zbuffer;
}

inline void plotter::drawLine(const pixel<2>& p1, const pixel<2>& p2, Uint32 col)
{
	/*
	 * Tegner en linje mellem to pixels
	 */

	drawLine(AFRUND(p1.x), AFRUND(p1.y), AFRUND(p2.x), AFRUND(p2.y), col);

	return;
}

inline void plotter::drawLineZ(const punkt& p1, const punkt& p2, Uint32 col)
{
	skalar len = (p2.pixl - p1.pixl).norm();	// Længden af linjen i pixels
	vektor<3> vec = p2.proj - p1.proj;			// Vektoren fra p1 til p2

	for (skalar t = 0; t <= 1; t += 1/len)
	{
		punkt r;
		r.proj = p1.proj + vec*t;	// Interpolér (også z-koordinaten)
		r.pixl = coord<2>(r.proj.x, r.proj.y).toPixel(this);

		int x = AFRUND(r.pixl.x);
		int y = AFRUND(r.pixl.y);

		if (!(0 < x && x < this->w)) continue;	// Pixel udenfor plot, skip den
		if (!(0 < y && y < this->h)) continue;	// ..

		// Tjek om der allerede er tegnet noget foran
		if (r.proj.z >= zbuffer->data[x][y]) {
			zbuffer->data[x][y] = r.proj.z;		// Opdatér Z-buffer
			AG_PUT_PIXEL2(this->plot, x, y, col);
		}
	}

	return;
}

inline Uint32 plotter::hsv(skalar h, skalar s, skalar v)
{
	/*
	 * Kilde:
	 * http://en.wikipedia.org/w/index.php?title=HSL_and_HSV&oldid=282264028#Conversion_from_HSV_to_RGB
	 */

	struct colour {
		skalar r, g, b;
		colour() {r = g = b = 1.0;} // Default til hvid
	} c;

	skalar f = h / 60;
	int hi = int(floor(f)) % 6;
	f -= floor(f);

	skalar p = v*(1 - s);
	skalar q = v*(1 - f*s);
	skalar t = v*(1 - (1-f)*s);

	switch(hi) {
		case 0:	c.r = v; c.g = t; c.b = p;	break;
		case 1:	c.r = q; c.g = v; c.b = p;	break;
		case 2:	c.r = p; c.g = v; c.b = t;	break;
		case 3:	c.r = p; c.g = q; c.b = v;	break;
		case 4:	c.r = t; c.g = p; c.b = v;	break;
		case 5:	c.r = v; c.g = p; c.b = q;	break;
	}

	c.r *= 255; // Gå fra [0; 1] til [0; 255]
	c.g *= 255; // ...
	c.b *= 255;

	// Pack farve komponenterne sammen til én værdi
	return (Uint32(c.b) << 16) | (Uint32(c.g) << 8) | Uint32(c.r);
}

inline Uint32 plotter::GetColour(skalar valdB, skalar mindB, skalar maxdB)
{
	/*
	 * Returnér farve svarende til val, lineært interpoleret mellem violet
	 * (minimum) og rød (maksimum).
	 * Det antages at val, min og max er angivet i decibel.
	 *
	 * For val > max går farven asymptotisk fra rød mod hvid for val --> Inf
	 * For val < min går farven asymptotisk fra violet mod sort for val --> -Inf
	 *
	 * NB: For nemmere at skelne farver som tilhører min-max intervallet, og dem
	 * som ligger udenfor, er overgangen gjort en smule diskontinuert.
	 */

	// Håndtér for singulariteter, Inf og 0 (= -Inf dB) feltstyrker
	if (!isfinite(valdB))
	{	/* For mere info jvf. http://www.johndcook.com/IEEE_exceptions_in_cpp.html */
	//	printf("_fpclass(valdB) = %i\n", _fpclass(valdB));
		if (isnan(valdB)) {
			return colorWhite;	// Singularitet (+Inf feltstyrke) = _FPCLASS_QNAN
		} else if(isinf(valdB)) {
			return colorBlack;	// 0 feltstyrke (-Inf dB) = _FPCLASS_NINF
		} else {
			return colorGrey;	// Ukendt
		}
	}

	/*
	 * Håndtér under/overflow ift. farveskalaens [min; max].
	 *
	 * Feltstyrkerne er angivet i decibel i henhold til 20*log10(|X|).
	 * Feltstyrken for vektorfelter spænder [0;Inf[, og er ikke negativ pga vektornorm.
	 * Feltstyrken for skalarpotentialet kan sagtens være negativ.
	 * Skalære felters styrke kan derfor ikke angives i en dB skala uden tab af fortegn!
	 */
	if (valdB > maxdB) {
	//	max = pow(10, maxdB/20.0);	// Gør lineær
	//	val = pow(10, valdB/20.0);	// ...
	//	skalar s = max/val;			// HSV saturation: 0.0 -> hvid; 1.0 -> rød
		skalar s = pow(10, (maxdB-valdB)/20); // Ækvivalent med ovenstående

		assert(0 <= s && s <=1);
		return hsv(0, s, 0.9);		// HSV hue = 0 -> rød
	} else if (valdB < mindB) {
	//	min = pow(10, mindB/20.0);	// Gør lineær
	//	val = pow(10, valdB/20.0);	// ...
	//	skalar v = val/min;			// HSV value: 0.0 -> sort; 1.0 -> violet
		skalar v = pow(10, (valdB-mindB)/20); // Ækvivalent med ovenstående

		assert(0 <= v && v <=1);
		return hsv(300, 0.9, v);	// HSV hue = 300 -> violet
	}

	// Interpolér farve lineært mellem decibel værdier
	skalar h = (valdB - mindB)/(maxdB - mindB);
	return hsv((1-h)*300, 1.0, 1.0);
}

inline void plotter::drawLine(int x1, int y1, int x2, int y2, Uint32 col)
{
	drawLine(this->plot, x1, y1, x2, y2, col); // default til at tegne i plottet

	return;
}

inline void plotter::drawLine(AG_Surface *surf, int x1, int y1, int x2, int y2, Uint32 col)
{
	/*
	 * Kilde:
	 * http://roguebasin.roguelikedevelopment.org/index.php?title=Bresenham%27s_Line_Algorithm&oldid=12887
	 */

	int delta_x = abs(x2 - x1) << 1;
	int delta_y = abs(y2 - y1) << 1;

	// if x1 == x2 or y1 == y2, then it does not matter what we set here
	signed char ix = x2 > x1 ? 1 : -1;
	signed char iy = y2 > y1 ? 1 : -1;

	AG_PUT_PIXEL2_CLIPPED(surf, x1, y1, col);

	if (delta_x >= delta_y) {
		// error may go below zero
		int error = delta_y - (delta_x >> 1);

		while (x1 != x2) {
			if (error >= 0) {
				if (error || (ix > 0)) {
					y1 += iy;
					error -= delta_x;
				}
			}

			x1 += ix;
			error += delta_y;

			AG_PUT_PIXEL2_CLIPPED(surf, x1, y1, col);
		}
	} else {
		// error may go below zero
		int error = delta_x - (delta_y >> 1);

		while (y1 != y2) {
			if (error >= 0) {
				if (error || (iy > 0)) {
					x1 += ix;
					error -= delta_y;
				}
			}

			y1 += iy;
			error += delta_x;

			AG_PUT_PIXEL2_CLIPPED(surf, x1, y1, col);
		}
	}

	return;
}

inline void plotter::pseudoLine(vector< list<int> >& ytable, const pixel<2>& p1, const pixel<2>& p2)
{
	int x1, y1, x2, y2;
	x1 = AFRUND(p1.x); y1 = AFRUND(p1.y);
	x2 = AFRUND(p2.x); y2 = AFRUND(p2.y);

	/*
	 * Kilde:
	 * http://roguebasin.roguelikedevelopment.org/index.php?title=Bresenham%27s_Line_Algorithm&oldid=12887
	 */

	int delta_x = abs(x2 - x1) << 1;
	int delta_y = abs(y2 - y1) << 1;

	// if x1 == x2 or y1 == y2, then it does not matter what we set here
	signed char ix = x2 > x1 ? 1 : -1;
	signed char iy = y2 > y1 ? 1 : -1;

	ytable[y1].push_back(x1);

	if (delta_x >= delta_y) {
		// error may go below zero
		int error = delta_y - (delta_x >> 1);

		while (x1 != x2) {
			if (error >= 0) {
				if (error || (ix > 0)) {
					y1 += iy;
					error -= delta_x;
				}
			}

			x1 += ix;
			error += delta_y;

			ytable[y1].push_back(x1);
		}
	} else {
		// error may go below zero
		int error = delta_x - (delta_y >> 1);

		while (y1 != y2) {
			if (error >= 0) {
				if (error || (iy > 0)) {
					x1 += ix;
					error -= delta_y;
				}
			}

			y1 += iy;
			error += delta_x;

			ytable[y1].push_back(x1);
		}
	}

	return;
}

inline void plotter::drawInOutVector(int x, int y, bool out, const Uint32 c)
{
	// Tegn en lille cirkel (5x5 pixel kvadrat med "afrundede" hjørner)
	for (int k = -1; k <= 1; k++) {
		AG_PUT_PIXEL2_CLIPPED(plot, x+k, y+2, c);
		AG_PUT_PIXEL2_CLIPPED(plot, x+k, y-2, c);
		AG_PUT_PIXEL2_CLIPPED(plot, x+2, y+k, c);
		AG_PUT_PIXEL2_CLIPPED(plot, x-2, y+k, c);
	}

	if (out) {
		AG_PUT_PIXEL2_CLIPPED(plot, x, y, c); // Vektor peger ud af skærmen
	} else {
		for (int k = -1; k <= 1; k++) {
			AG_PUT_PIXEL2_CLIPPED(plot, x+k, y, c);
			AG_PUT_PIXEL2_CLIPPED(plot, x, y+k, c);
		}
	}

	AG_PUT_PIXEL2_CLIPPED(plot, x, y, colorWhite);

	return;
}

inline bool plotter::drawNormedVector(const coord<2>& start, const vektor<2>& vek, Uint32 col)
{
	/*
	 * Plotter en norm-standardiseret vektor i koordinatsystemet.
	 * Koordinaterne i start og vek, er i alm koordinater (ikke pixels).
	 */

	vektor<2> vekny = (vek / vek.norm()) * start.norm(); // Undgå trunkeringsfejl
	coord<2> slut = start + vekny;			// Coord for vek's hoved

	pixel<2> pStart = start.toPixel(this);	// Pixel hvor vek begynder
	pixel<2> pSlut = slut.toPixel(this);	// Pixel hvor vek slutter
	vektor<2> vPix(pStart, pSlut);			// Pixel-vektor ækvivalenten til vek

	vPix = vPix / vPix.norm();				// Normér ækvivalenten til 1 pixel
	vPix = vPix * pilLen;					// Pilen skalerer ikke med plot size

	pixel<2> pSlut2 = pStart + vPix;		// Pixel hvor normeret vek slutter
	return drawArrow(pStart, pSlut2, col);
}

inline void plotter::drawNormedVectorZ(	const punkt& start, const vektor<3>& vek,
										skalar vNorm, Uint32 col)
{
	vektor<3> vUnit = vek / vek.norm();
	vektor<3> v0 = vUnit * vNorm;

	coord<3> cslut = start.orig + v0;
	punkt slut = cslut.toPunkt(this);

	drawLineZ(start, slut, col);

	return;
}

inline bool plotter::drawArrow(const pixel<2>& pStart, const pixel<2>& pSlut, Uint32 col)
{
	/*
	 * Tegner en pil fra pStart til pSlut, hvor disse er angivet i pixels.
	 * Pilens hoved tegnes ved pSlut.
	 */

	// Vi gider ikke at tegne pile som ikke er fuldkommen indeholdt i plottet
	if (!(0 <= pStart.x && pStart.x <= w))	return false; // pStart skal være inde i plottet
	if (!(0 <= pStart.y && pStart.y <= h)) 	return false;
	if (!(0 <= pSlut.x && pSlut.x <= w)) 	return false; // pSlut skal være inde i plottet
	if (!(0 <= pSlut.y && pSlut.y <= h)) 	return false;

	vektor<2> v0(pStart, pSlut);		// Pixel-vektoren fra pstart til pSlut
	vektor<2> vn = v0 / v0.norm();		// Normaliseret v0: 1 pixel lang
	vektor<2> vt = vn.orto()*headWid;	// Vektor ud til hovedets hjørner

	// Hovedet udspændes af {pA, pB, pSlut}
	pixel<2> pH = pSlut - vn*headLen;	// Midtpunkt for hovedets bagstykke
	pixel<2> pA = pH + vt;				// Pixel-koordinat til hjørne A
	pixel<2> pB = pH - vt;				// Pixel-koordinat til hjørne B

	drawLine(pStart, pH, col);	// Pilens "krop"; kan også gå fra pStart til pSlut
	drawLine(pA, pSlut, col);	// Venstre/højre linjestykke af hovedet
	drawLine(pB, pSlut, col);	// Højre/venstre linjestykke af hovedet
	drawLine(pB, pA, col);		// Bagstykket af hovedet

	// Markér pixel hvor feltet er målt i (= start af vektor-pilen i det punkt)
	AG_PUT_PIXEL2_CLIPPED(plot, AFRUND(pStart.x), AFRUND(pStart.y), colorWhite);

	return true;	// Succes, pilen blev tegnet
}

inline void plotter::drawAxis(char axis1, char axis2)
{
	/*
	 *	Tegner akser i 2D (og 3D).
	 */

	AG_Surface *tmp;
	pixel<2> pxBeg = coord<2>(xmin, 0).toPixel(this);
	pixel<2> pxEnd = coord<2>(xmax, 0).toPixel(this);
	pixel<2> pyBeg = coord<2>(0, ymin).toPixel(this);
	pixel<2> pyEnd = coord<2>(0, ymax).toPixel(this);

	pxBeg.x = 1; pxEnd.x = w-1;	// Fiksér akserne (undgå jitter @ AFRUND)
	pyBeg.y = h-1; pyEnd.y = 1;

	// Plot 2D akser
	drawArrow(pxBeg, pxEnd);	// Tegn 1. aksen
	drawArrow(pyBeg, pyEnd);	// Tegn 2. aksen

	// Afstanden mellem tickmarks er altid en 10'er potens; og lavere end max range
	skalar unit = ((xrange>yrange) ? xrange : yrange); // Find max range,
	unit = pow(10, POTENS10(unit) - 1);	// rund ned til nærmeste 10'er potens,
										// -1: og lav indelingen en tak mere fin

	// Tegn tickmarks
	for (skalar i = 0; i < xmax; i += unit) {
		pixel<2> tic = coord<2>(i, 0).toPixel(this);
		if (tic.x > w-headLen-1) break; // Tegn ikke tickmarks inde i pilens hoved
		drawLine(AFRUND(tic.x), AFRUND(tic.y)+tickLen,
				 AFRUND(tic.x), AFRUND(tic.y)-tickLen);
	}
	for (skalar i = 0; i > xmin; i -= unit) {
		pixel<2> toc = coord<2>(i, 0).toPixel(this);
		drawLine(AFRUND(toc.x), AFRUND(toc.y)+tickLen,
				 AFRUND(toc.x), AFRUND(toc.y)-tickLen);
	}
	for (skalar i = 0; i < ymax; i += unit) {
		pixel<2> tic = coord<2>(0, i).toPixel(this);
        if (tic.y < headLen+1) break; // Tegn ikke tickmarks inde i pilens hoved
		drawLine(AFRUND(tic.x)+tickLen, AFRUND(tic.y),
				 AFRUND(tic.x)-tickLen, AFRUND(tic.y));
	}
	for (skalar i = 0; i > ymin; i -= unit) {
		pixel<2> toc = coord<2>(0, i).toPixel(this);
		drawLine(AFRUND(toc.x)+tickLen, AFRUND(toc.y),
				 AFRUND(toc.x)-tickLen, AFRUND(toc.y));
	}

	// Tegn labels
	tmp = AG_TextRenderf("%c-akse", axis1);
	AG_SurfaceBlit(tmp, NULL, this->plot,
		AFRUND(pxEnd.x - tmp->w),
		AFRUND(pxEnd.y - tmp->h - headWid));
	{
		skalar potens;

		potens = POTENS10(xmax);
		tmp = LABEL(xmax, potens);
		AG_SurfaceBlit(tmp, NULL, this->plot,
			AFRUND(pxEnd.x - tmp->w),
			AFRUND(pxEnd.y + headWid));

		potens = POTENS10(xmin);
		tmp = LABEL(xmin, potens);
		AG_SurfaceBlit(tmp, NULL, this->plot,
			AFRUND(pxBeg.x),
			AFRUND(pxBeg.y + headWid));
	}

	tmp = AG_TextRenderf("%c-akse", axis2);
	AG_SurfaceBlit(tmp, NULL, this->plot,
		AFRUND(pyEnd.x + headLen),
		AFRUND(pyEnd.y));
	{
		skalar potens;

		potens = POTENS10(ymax);
		tmp = LABEL(ymax, potens);
		AG_SurfaceBlit(tmp, NULL, this->plot,
			AFRUND(pyEnd.x - tmp->w - headWid),
			AFRUND(pyEnd.y));

		potens = POTENS10(ymin);
		tmp = LABEL(ymin, potens);
		AG_SurfaceBlit(tmp, NULL, this->plot,
			AFRUND(pyBeg.x - tmp->w - headWid),
			AFRUND(pyBeg.y - tmp->h));
	}

	return;
}

void plotter::msgNoSuitableObjects()
{
	AG_TextColor32(colorWhite);
	AG_Surface *tmp = AG_TextRenderf(	"The selected field is not applicable\n"
										"to any of the currently added objects.\n"
										"\n"
										"You can add suitable objects via the\n"
										"\"Add/remove objects\" menu.");
	AG_SurfaceBlit(tmp, NULL, plot,
		AFRUND((plot->w - tmp->w)/2),
		AFRUND((plot->h - tmp->h)/2));

	return;
}

void plotter::msgNoObjects()
{
	AG_TextColor32(colorWhite);
	AG_Surface *tmp = AG_TextRenderf(	"No objects currently added.\n"
										"\n"
										"You can add objects via the\n"
										"\"Add/remove objects\" menu.");
	AG_SurfaceBlit(tmp, NULL, plot,
		AFRUND((plot->w - tmp->w)/2),
		AFRUND((plot->h - tmp->h)/2));

	return;
}

inline void plotter::calcPlaneZ(const Medium& vem, const vector<punkt>& punkter,
								map< int,map<int,skalar> >& intensities)
{
	int ymin, ymax;
	ymin = ymax = AFRUND(punkter[0].pixl.y);

	for (int k = 0; k < punkter.size(); k++)
	{
		int tmp = AFRUND(punkter[k].pixl.y);

		if (ymin > tmp) ymin = tmp;
		if (ymax < tmp) ymax = tmp;
	}

	vector< list<int> > ytable;
	ytable.resize(this->h, list<int>());

	assert(punkter.size() == 4);
	pseudoLine(ytable, punkter[0].pixl, punkter[1].pixl);
	pseudoLine(ytable, punkter[1].pixl, punkter[2].pixl);
	pseudoLine(ytable, punkter[2].pixl, punkter[3].pixl);
	pseudoLine(ytable, punkter[3].pixl, punkter[0].pixl);

	vektor<3> v1(punkter[0].offrot, punkter[1].offrot);
	vektor<3> v2(punkter[0].offrot, punkter[3].offrot);
	vektor<3> normal = v1.kryds(v2); // Normalvektoren til fladen

	for (int y = ymin; y <= ymax; y++)
	{
		ytable[y].sort(); // Højst 3 elementer heri, så det går nok...

		int x1 = ytable[y].front();
		int x2 = ytable[y].back();

		map<int, skalar>& ycurrent = intensities[y]; // Cache map-lookup

		for (int x = x1; x <= x2; x++)
		{
			if (!(0 < x && x < this->w)) continue; // Pixel udenfor plot, skip den
			if (!(0 < y && y < this->h)) continue; // ...

			// Find tilhørende 2D projektionskoordinater for vores pixel
			coord<2> pm = pixel<2>(x, y).toCoord(this);

			// Find 3D punktet r.offrot som, når projiceret, giver pm
			punkt r; // Vores feltsamplingspunkt
			skalar G = f - W;
			vektor<3> S0((pm.x*f)/G, (pm.y*f)/G, 0);
			vektor<3> u(-pm.x/G, -pm.y/G, 1);
			skalar t = normal.dot(punkter[0].offrot - S0)/normal.dot(u);
			r.offrot = S0 + u*t;

			// Afgør om der allerede er tegnet noget foran vores punkt
			skalar zdepth = r.offrot.z / (f - r.offrot.z);
			if ( !(zdepth > zbuffer->data[x][y]) )
				continue;	// Der er noget foran, så vi gider slet ikke at regne
							// feltet ud når brugeren alligevel ikke kan se det

			// Der er frit udsyn; opdatér z-buffer (og regn feltet ud, plot det)
			zbuffer->data[x][y] = zdepth;

			// Rotér og translatér tilbage til korrekt samplingspunkt
			r.orig = r.offrot.rotz(-alpha);	// Rotér tilbage i omvendt rækkefølge
			r.orig = r.orig.rotx(-beta);	// ...
			r.orig = r.orig.rotz(-gamma);
			r.orig = r.orig + centrum;		// Læg offset til

			// Beregn resulterende feltvektor i samplingspunktet
			fieldPoint totalField;
			for (int i = 0; i < vem.size(); i++) {
				totalField = totalField + vem[i]->felt[vem.mainField]->calc(r.orig);
			}

			// Gem intensiteten af beregnet feltvektor
			ycurrent[x] = 20*log10(totalField.getItensity()); // = intensities[y][x] = dB
		}
	}

	return;
}

void plotter::tegn3d(const Medium& vem, skalar f, skalar W, coord<3> centrum)
{
	/*
	 * Tegner en kube som repræsenterer det 3D rum hvori feltet ønskes undersøgt,
	 * samt 3 flader hvorpå den totale feltintensitet er farvekodet.
	 */

	map< int,map<int,skalar> > intensities;
	zbuffer = new matrix<>(this->w, this->h, -INF);

	this->f = f;
	this->W = W;
	this->centrum = centrum;

	this->alpha = param.view3d.alpha * (pi/180) + EPS;	// Konvertér til radianer
	this->beta  = param.view3d.beta * (pi/180) + EPS;	// EPS: Undgå floatingpoint quirks
	this->gamma = param.view3d.gamma * (pi/180) + EPS;	// ...


	polygon box;

	viewbox3d tmp = param.view3d;
	box.addPoint("MIN",  coord<3>(tmp.x.min, tmp.y.min, tmp.z.min).toPunkt(this)); // Bund
	box.addPoint("zmax", coord<3>(tmp.x.min, tmp.y.min, tmp.z.max).toPunkt(this)); // Bund
	box.addPoint("ymax", coord<3>(tmp.x.min, tmp.y.max, tmp.z.min).toPunkt(this)); // Top
	box.addPoint("xmin", coord<3>(tmp.x.min, tmp.y.max, tmp.z.max).toPunkt(this)); // Top
	box.addPoint("xmax", coord<3>(tmp.x.max, tmp.y.min, tmp.z.min).toPunkt(this)); // Bund
	box.addPoint("ymin", coord<3>(tmp.x.max, tmp.y.min, tmp.z.max).toPunkt(this)); // Bund
	box.addPoint("zmin", coord<3>(tmp.x.max, tmp.y.max, tmp.z.min).toPunkt(this)); // Top
	box.addPoint("MAX",  coord<3>(tmp.x.max, tmp.y.max, tmp.z.max).toPunkt(this)); // Top

	Uint32 colorAxis = colorGrey;
	box.makeLine("MIN",  "zmax", colorAxis);	// Bund
	box.makeLine("zmax", "ymin", colorAxis);	// ...
	box.makeLine("ymin", "xmax", colorAxis);
	box.makeLine("xmax", "MIN",  colorAxis);
	box.makeLine("MAX",  "xmin", colorAxis);	// Top
	box.makeLine("xmin", "ymax", colorAxis);	// ...
	box.makeLine("ymax", "zmin", colorAxis);
	box.makeLine("zmin", "MAX",  colorAxis);
	box.makeLine("zmax", "xmin", colorAxis);	// Ventre side
	box.makeLine("ymax", "MIN",  colorAxis);	// ...
	box.makeLine("xmax", "zmin", colorAxis);	// Højre side
	box.makeLine("MAX",  "ymin", colorAxis);	// ...

	box.makeFlade(this, "MIN", "zmax", "xmin", "ymax", 0, slice.planes[0]); // x-aksen er normal herpå
	box.makeFlade(this, "MIN", "zmax", "ymin", "xmax", 1, slice.planes[1]); // y-aksen ...
	box.makeFlade(this, "MIN", "ymax", "zmin", "xmax", 2, slice.planes[2]); // z-aksen ...

	box.addPoint("origo", coord<3>(0, 0, 0).toPunkt(this));
	box.addPoint("xunit", coord<3>(1, 0, 0).toPunkt(this));
	box.makeLine("origo", "xunit", colorRed);	// unit vector langs x-aksen
	box.addPoint("yunit", coord<3>(0, 1, 0).toPunkt(this));
	box.makeLine("origo", "yunit", colorGreen);	// unit vector langs y-aksen
	box.addPoint("zunit", coord<3>(0, 0, 1).toPunkt(this));
	box.makeLine("origo", "zunit", colorBlue);	// unit vector langs z-aksen

	for (int i = 0; i < box.linier.size(); i++) {
		punkt p1 = *box.linier[i].par.first;
		punkt p2 = *box.linier[i].par.second;

		drawLineZ(p1, p2, box.linier[i].farve); // Plot linjer
	}

	for (int i = 0; i < box.flader.size(); i++)	// Beregn feltintensitet
		calcPlaneZ(vem, box.flader[i].punkter, intensities);
	if (param.adaptiveColor == true)
		adaptivFarve3D(&intensities);	// Find passende range for farveskalaen
	plotBackground3D(&intensities);		// Vis feltintensiteten på fladerne


	// Tegn farveskalaen
	for (int i = 0; i < skala->w; i++)
		drawLine(skala, i, 0, i,skala->h-1, GetColour(i, 0, skala->w));

	return;
}

void plotter::calcIntensity(const Medium& vem, matrix<>& m, Field::pripot felt)
{
	/*
	 * Beregner feltintensiteten af det totale vektorfelt (ikke projiceret)
	 */

	int A, B, C;
	switch(param.cslice) {
		case 'x': A = 0; B = 1; C = 2; break;
		case 'y': A = 1; B = 2; C = 0; break;
		case 'z': A = 2; B = 0; C = 1; break;
	}

	// Loop alle pixels i gennem, find intensiteten af totalfeltet
	for (int y = 1; y < h; y++) {
	for (int x = 1; x < w; x++) {
			fieldPoint totalField;
			coord<2> c2d = pixel<2>(x, y).toCoord(this);
			coord<3> c3d;

			// Projicer 2D skærm-punkt til ægte 3D punkt hvori feltet skal beregnes
			c3d.data[A][0] = param.nslice;
			c3d.data[B][0] = c2d.x;
			c3d.data[C][0] = c2d.y;

			for (int i = 0; i < vem.size(); i++) {
				totalField = totalField + vem[i]->felt[felt]->calc(c3d);
			}

			m.data[x][y] = 20*log10(totalField.getItensity());
		}
	}

	return;
}

void plotter::adaptivFarve(const matrix<>& m)
{
	vector<skalar> liste;
	liste.reserve(w*h);

	for (int x = 1; x <= w-1; x++) {
	for (int y = 1; y <= h-1; y++) {
			skalar val = m.data[x][y];
			if (!isfinite(val))		// Spring Inf's og NaN's over
				continue;
			liste.push_back(val);
		}
	}

	if (liste.empty()) {
		// Dette er tilfældet hvis et objekt er tilføjet, men har styrke lig 0.
		// Med styrke lig 0, bliver feltintensiteten også 0 = -Inf dB.
		farve.max = -INF;
		farve.min = -INF;
		return;
	}

	sort(liste.begin(), liste.end()); // Sortér efter stigende værdi

	int nth = int(liste.size()*0.99);
	int mth = int(liste.size()*0.05);
	assert(0 <= nth && nth < liste.size());
	assert(0 <= mth && mth < liste.size());
	farve.max = liste[nth];		// Find 99% største værdi
	farve.min = liste[mth];		// Find  5% største værdi

	return;
}

void plotter::adaptivFarve3D(map< int,map<int,skalar> > *intensities)
{
	map< int,map<int,skalar> >::iterator ity;
	map< int,skalar>::iterator itx;
	vector<skalar> liste;
//	liste.reserve();

	for (ity = intensities->begin(); ity != intensities->end(); ity++) {
	for (itx = ity->second.begin(); itx != ity->second.end(); itx++) {
			skalar val = itx->second;

			if (!isfinite(val))		// Spring Inf's og NaN's over
				continue;
			liste.push_back(val);
		}
	}

	if (liste.empty()) {
		// Dette er tilfældet hvis et objekt er tilføjet, men har styrke lig 0.
		// Med styrke lig 0, bliver feltintensiteten også 0 = -Inf dB.
		farve.max = -INF;
		farve.min = -INF;
		return;
	}

	sort(liste.begin(), liste.end()); // Sortér efter stigende værdi

	int nth = int(liste.size()*0.99);
	int mth = int(liste.size()*0.05);
	assert(0 <= nth && nth < liste.size());
	assert(0 <= mth && mth < liste.size());
	farve.max = liste[nth];		// Find 99% største værdi
	farve.min = liste[mth];		// Find  5% største værdi

	return;
}

void plotter::plotBackground(const matrix<>& m)
{
	// Plot feltintensiteten af det totale vektorfelt (ikke projiceret) som baggrund
	for (int y = 1; y < h; y++) {
	for (int x = 1; x < w; x++)	{
			AG_PUT_PIXEL2(plot, x, y, GetColour(m.data[x][y], farve.min, farve.max));
		}
	}

	return;
}

void plotter::plotBackground3D(map< int,map<int,skalar> > *intensities)
{
	map< int,map<int,skalar> >::iterator ity;
	map< int,skalar>::iterator itx;

	for (ity = intensities->begin(); ity != intensities->end(); ity++) {
		int y = ity->first;

		for (itx = ity->second.begin(); itx != ity->second.end(); itx++) {
			int x = itx->first;
			skalar val = itx->second;

			Uint32 col = GetColour(val, farve.min, farve.max);
			AG_PUT_PIXEL2(plot, x, y, col);
		}
	}

	return;
}

void plotter::plotEquipotentialLines(const Medium& vem)
{
	matrix<> mm(w, h);			// Gemmer feltintensiteten for hver pixel

	calcIntensity(vem, mm, Field::potential); // Beregn V-felt for hver eneste pixel

	// Find intervaller
	vector<skalar> liste;
	liste.reserve(w*h);

	for (int x = 1; x <= w-1; x++) {
	for (int y = 1; y <= h-1; y++) {
			skalar val = mm.data[x][y];
			if (!isfinite(val))		// Spring Inf's og NaN's over
				continue;
			liste.push_back(val);
		}
	}

	sort(liste.begin(), liste.end()); // Sortér efter stigende værdi

	const float inddel[] = {0.20, 0.40, 0.60, 0.80};
	const int ndel = 4;
	skalar step[ndel];

	for (int i = 0; i < ndel; i++) {
		step[i] = liste[int(liste.size()*inddel[i])];
	}

	for (int y = 1; y < h; y++) {
	for (int x = 1; x < w; x++)	{
			skalar val = mm.data[x][y];

			int t = ndel;
			for (int k = 0; k < ndel; k++) {
				if (val < step[k]) {
					t = k;
					break;
				}
			}

			mm.data[x][y] = t;
		}
	}

	for (int y = 2; y < h-1; y++) {
	for (int x = 2; x < w-1; x++) {
			skalar delta1 = int(mm.data[x][y]) - int(mm.data[x+1][y]);
		    skalar delta2 = int(mm.data[x][y]) - int(mm.data[x][y+1]);

			if (delta1 < 0)
				AG_PUT_PIXEL2(plot, x+1, y, colorWhite);
			if (delta2 < 0)
				AG_PUT_PIXEL2(plot, x, y+1, colorWhite);
			if (delta1 > 0 || delta2 > 0)
				AG_PUT_PIXEL2(plot, x, y, colorWhite);
		}
	}

	return;
}

void plotter::plotVectorField(const Medium& vem)
{
	int A, B, C;
	switch(param.cslice) {
		case 'x': A = 0; B = 1; C = 2; break;
		case 'y': A = 1; B = 2; C = 0; break;
		case 'z': A = 2; B = 0; C = 1; break;
	}

	// Plot det totale vektorfeltet projiceret ind på slice-planet
	for (int y = 2; y < h; y += distPile) {
	for (int x = 2; x < w; x += distPile) {
			fieldPoint totalField;
			coord<2> c2d = pixel<2>(x, y).toCoord(this); // Pixel til første/anden-akse koordinater
			coord<3> c3d;	// Koordinater til feltsampling i 3D-space (ligger i sliceplanet)
			vektor<2> v2d;	// Feltvektor projiceret ind på sliceplanet

			// Projicer 2D skærm-punkt til det ægte 3D punkt hvori feltet skal beregnes
			c3d.data[A][0] = param.nslice;
			c3d.data[B][0] = c2d.x;
			c3d.data[C][0] = c2d.y;

			// Beregn resulterende feltvektor i 3D punktet
			for (int i = 0; i < vem.size(); i++) {
				totalField = totalField + vem[i]->felt[vem.mainField]->calc(c3d);
			}

			// Og projicer feltvektor ind på skærmplanet
			v2d.x = totalField.vfelt.data[B][0];
			v2d.y = totalField.vfelt.data[C][0];

			// Vis projiceret feltvektor på skærmen
			drawNormedVector(c2d, v2d);
			if (v2d.norm() == 0)
				drawInOutVector(x, y, totalField.vfelt.data[A][0] > 0);
		}
	}

	return;
}

void plotter::tegn2d(const Medium& vem)
{
	/*
	 * Plotter givne em-objekter, hvis de kan; ellers skriv fejl medd i selve plottet.
	 */

	matrix<> m(w, h);			// Gemmer feltintensiteten for hver pixel

	calcIntensity(vem, m, vem.mainField); // Beregn feltintensiteten for hver eneste pixel
	if (param.adaptiveColor == true)
		adaptivFarve(m);		// Find passende range for farveskalaen
	plotBackground(m);			// Vis feltintensiteten som baggrund

	if (vem[0]->felt[!vem.mainField]->type == fieldPoint::skalar_felt)
		plotEquipotentialLines(vem);	// Plot udvalgte ækvipotentiallinjer
	if (vem[0]->felt[vem.mainField]->type == fieldPoint::vektor_felt)
		plotVectorField(vem);			// Plot vektorfeltet

	// Akserne er altid synlige, så de tegnes sidst
	switch(param.cslice) {
		case 'x': drawAxis('y', 'z'); break;
		case 'y': drawAxis('z', 'x'); break;
		case 'z': drawAxis('x', 'y'); break;
	}


	// Tegn farveskalaen
	for (int i = 0; i < skala->w; i++)
		drawLine(skala, i, 0, i,skala->h-1, GetColour(i, 0, skala->w));

	return;
}

AG_Surface* plotter::getIntensity()
{
	return this->skala;
}

AG_Surface* plotter::getPlot()
{
	return this->plot;
}
